The example explains effective permeabilty calculations using PoreSpy and OpenPNM software. The simulation is performed on X-ray tomography image of BereaSandstone. The calculated effective permeablity value can compared with value report in Dong et al.
In [1]:
import os
import imageio
import scipy as sp
import numpy as np
import openpnm as op
import porespy as ps
import matplotlib.pyplot as plt
np.set_printoptions(precision=4)
np.random.seed(10)
%matplotlib inline
Give path to image file and load the image. Please note image should be binarized or in boolean format before performing next steps.
In [2]:
path = '../fixtures/ICL-Sandstone(Berea)/'
file_format = '.tif'
file_name = 'Berea'
file = file_name + file_format
fetch_file = os.path.join(path, file)
im = imageio.mimread(fetch_file)
im = ~np.array(im, dtype=bool)
In [3]:
#NBVAL_IGNORE_OUTPUT
fig, ax = plt.subplots(1, 3, figsize=(12,5))
ax[0].imshow(im[:, :, 100]);
ax[1].imshow(ps.visualization.show_3D(im));
ax[2].imshow(ps.visualization.sem(im));
ax[0].set_title("Slice No. 100 View");
ax[1].set_title("3D Sketch");
ax[2].set_title("SEM View");
In [4]:
print(ps.metrics.porosity(im))
The SNOW algorithm (an accronym for Sub-Network from an Oversegmented Watershed) was presented by Gostick. The algorithm was used to extract pore network from BereaSandstone image.
In [5]:
#NBVAL_IGNORE_OUTPUT
resolution = 5.345e-6
net = ps.networks.snow(im=im, voxel_size=resolution)
The output from the SNOW algorithm above is a plain python dictionary containing all the extracted pore-scale data, but it is NOT yet an OpenPNM network. We need to create an empty network in OpenPNM, then populate it with the data from SNOW:
In [6]:
ws = op.Workspace()
proj = op.Project()
pn = op.network.GenericNetwork(name=file_name, project=proj)
pn.update(net) # Fills 'pn' with data from 'net'
Now we can print the network to see how the transferred worked:
In [7]:
print(pn)
Remove isolated pores or cluster of pores from the network by checking it network health. Make sure ALL keys in network health functions have no value.
In [8]:
h = pn.check_network_health()
op.topotools.trim(network=pn, pores=h['trim_pores'])
h = pn.check_network_health()
print(h)
Assign conduit shape to calculate flow from pore i to pore j. In this simulation spherical pore and cylinderical throats geometry is considered.
In [9]:
mod = op.models.geometry.throat_endpoints.spherical_pores
pn.add_model(propname='throat.endpoints', model=mod)
mod = op.models.geometry.throat_length.conduit_lengths
pn.add_model(propname='throat.conduit_lengths', model=mod)
mod = op.models.geometry.pore_area.sphere
pn.add_model(propname='pore.area', model=mod)
In [10]:
geo = op.geometry.GenericGeometry(network=pn, pores=pn.Ps, throats=pn.Ts)
In this example air is considered as fluid passing through porous channels.
In [11]:
air = op.phases.Air(network=pn)
In [12]:
phys_air = op.physics.Standard(network=pn, phase=air, geometry=geo)
Select stokes flow algorithm for simulation and assign dirichlet boundary conditions in top and bottom faces of the network.
In [13]:
perm = op.algorithms.StokesFlow(network=pn, project=proj)
perm.setup(phase=air)
perm.set_value_BC(pores=pn.pores('top'), values=0)
perm.set_value_BC(pores=pn.pores('bottom'), values=101325)
perm.run()
air.update(perm.results())
Caclulate effective permeablity using hagen poiseuille equation. Use cross section area and flow length manually from image dimension.
In [14]:
resolution = 5.345e-6
Q = perm.rate(pores=pn.pores('bottom'), mode='group')
A = (im.shape[0] * im.shape[1]) * resolution**2
L = im.shape[2] * resolution
mu = air['pore.viscosity'].max()
delta_P = 101325 - 0
K = Q * L * mu / (A * delta_P)
print('The value of K is:', K/0.98e-12*1000, 'mD')